from __future__ import annotations
from dataclasses import dataclass
from typing import Tuple, List, Any, Callable, Optional
import numpy as np
import math

# --------------------------
# 实用小工具
# --------------------------
def _unit(v: np.ndarray) -> np.ndarray:
    v = np.asarray(v, float)
    n = np.linalg.norm(v)
    return v if n == 0 else v / n

def _merge_intervals(intervals: List[Tuple[float, float]], eps: float = 1e-10) -> List[Tuple[float, float]]:
    if not intervals:
        return []
    intervals = sorted(intervals)
    merged = [[intervals[0][0], intervals[0][1]]]
    for a, b in intervals[1:]:
        if a > merged[-1][1] + eps:
            merged.append([a, b])
        else:
            merged[-1][1] = max(merged[-1][1], b)
    return [(float(a), float(b)) for a, b in merged]

def _poly_cross_linear_linear_sqnorm(Wc, F, U0, D) -> np.ndarray:
    """
    || (Wc + F t) × (U0 + D t) ||^2 的四次多项式系数（降幂）[k4,k3,k2,k1,k0]
    """
    A0 = np.cross(Wc, U0)
    A1 = np.cross(Wc, D) + np.cross(F, U0)
    A2 = np.cross(F, D)
    k0 = float(np.dot(A0, A0))
    k1 = float(2.0 * np.dot(A0, A1))
    k2 = float(np.dot(A1, A1) + 2.0 * np.dot(A0, A2))
    k3 = float(2.0 * np.dot(A1, A2))
    k4 = float(np.dot(A2, A2))
    return np.array([k4, k3, k2, k1, k0], dtype=float)

def _poly_u2(U0, D) -> np.ndarray:
    """||U0 + D t||^2 的二次多项式系数（降幂）[c2,c1,c0]"""
    c2 = float(np.dot(D, D))
    c1 = float(2.0 * np.dot(U0, D))
    c0 = float(np.dot(U0, U0))
    return np.array([c2, c1, c0], dtype=float)

def _poly_dot_quadratic(Wc, F, U0, D) -> np.ndarray:
    """(Wc + F t)·(U0 + D t) 的二次多项式系数（降幂）[d2,d1,d0]"""
    d2 = float(np.dot(F, D))
    d1 = float(np.dot(Wc, D) + np.dot(F, U0))
    d0 = float(np.dot(Wc, U0))
    return np.array([d2, d1, d0], dtype=float)

def _poly_subtract_scaled_quadratic(quart, quad, scale) -> np.ndarray:
    """quart(t) - scale * quad(t) 的四次系数"""
    k4, k3, k2, k1, k0 = quart
    c2, c1, c0 = quad
    return np.array([k4, k3, k2 - scale * c2, k1 - scale * c1, k0 - scale * c0], dtype=float)

def _real_roots_in_interval(coeffs: np.ndarray, a: float, b: float, eps: float = 1e-12) -> list:
    """返回多项式（降幂）在 [a,b] 内的全部实根（排序）。自动裁剪前导零。"""
    c = np.array(coeffs, float)
    i = 0
    while i < len(c) - 1 and abs(c[i]) < eps:
        i += 1
    c = c[i:]
    if len(c) <= 1:
        return []
    roots = np.roots(c)
    real = []
    for r in roots:
        if abs(r.imag) < 1e-10:
            x = float(r.real)
            if a - 1e-12 <= x <= b + 1e-12:
                real.append(x)
    real.sort()
    return real

# --------------------------
# 解析评估（修正版，含端点）
# --------------------------
def intervals_point_los_analytic_fixed(
    M0: np.ndarray,
    aim_to: np.ndarray,        # 用于计算单位方向；若你已是单位向量也可直接传 dir
    Vm_speed: float,
    T: np.ndarray,             # 代表点：center/top/bottom 三选一
    C_det: np.ndarray,         # 起爆时刻的云心
    sink: float,
    t_det: float,
    r: float,
    t0: float,
    t1: float,
    eps: float = 1e-12
) -> tuple[List[Tuple[float, float]], float]:
    """
    代表点 LoS 的解析遮蔽区间。
    遮蔽条件等价于“线段-球相交”：(直线距离≤r 且 垂足在段内) 或 端点在球内。
    """
    if t1 <= t0:
        return [], 0.0

    # 线性模型
    M0 = np.asarray(M0, float)
    T  = np.asarray(T, float)
    Vm_dir = _unit(np.asarray(aim_to, float) - M0)  # ← 关键修正：传 aim_to，用此处单位化
    D = Vm_speed * Vm_dir
    U0 = M0 - T

    F = np.array([0.0, 0.0, -float(sink)], float)
    Wc = (np.asarray(C_det, float) - T) - F * float(t_det)  # 使得 W(t)=Wc + F t

    # 端点二次
    # Et(t) = ||W(t)||^2 - r^2
    Et = np.array([float(np.dot(F, F)),
                   float(2.0 * np.dot(Wc, F)),
                   float(np.dot(Wc, Wc) - r * r)], dtype=float)
    # Em(t) = ||C(t) - M(t)||^2 - r^2 = ||(C_det - F t_det - M0) + (F - D) t||^2 - r^2
    Qc = (np.asarray(C_det, float) - F * float(t_det)) - M0
    G  = F - D
    Em = np.array([float(np.dot(G, G)),
                   float(2.0 * np.dot(Qc, G)),
                   float(np.dot(Qc, Qc) - r * r)], dtype=float)

    # 直线距离 + 段内
    quart = _poly_cross_linear_linear_sqnorm(Wc, F, U0, D)
    u2    = _poly_u2(U0, D)
    g     = _poly_subtract_scaled_quadratic(quart, u2, r * r)       # 四次：||W×U||^2 - r^2||U||^2 ≤ 0
    h     = _poly_dot_quadratic(Wc, F, U0, D)                       # 二次：W·U
    h_le  = np.array([h[0] - u2[0], h[1] - u2[1], h[2] - u2[2]])    # 二次：W·U - ||U||^2 ≤ 0

    # 收集所有临界点
    crit = [t0, t1]
    crit += _real_roots_in_interval(g, t0, t1)
    crit += _real_roots_in_interval(h, t0, t1)
    crit += _real_roots_in_interval(h_le, t0, t1)
    crit += _real_roots_in_interval(Et, t0, t1)
    crit += _real_roots_in_interval(Em, t0, t1)
    crit = sorted(set(float(x) for x in crit))

    def _val_poly_desc(c: np.ndarray, t: float) -> float:
        v = 0.0
        for a in c:
            v = v * t + a
        return v

    def _cond(t: float) -> bool:
        # 直线 + 段内
        ok_line = (_val_poly_desc(g, t) <= eps) and (h[0]*t*t + h[1]*t + h[2] >= -eps) and (h_le[0]*t*t + h_le[1]*t + h_le[2] <= eps)
        # 端点命中
        ok_T = (Et[0]*t*t + Et[1]*t + Et[2]) <= eps
        ok_M = (Em[0]*t*t + Em[1]*t + Em[2]) <= eps
        return ok_line or ok_T or ok_M

    intervals: List[Tuple[float, float]] = []
    for i in range(len(crit) - 1):
        a, b = crit[i], crit[i+1]
        if b <= a:
            continue
        m = 0.5 * (a + b)
        if _cond(m):
            intervals.append((a, b))

    intervals = _merge_intervals(intervals)
    total = float(sum(b - a for a, b in intervals))
    return intervals, total

# --------------------------
# 可选：用于自检的小型扫描器（代表点 LoS）
# --------------------------
def segment_intersects_sphere(A: np.ndarray, B: np.ndarray, center: np.ndarray, r: float, eps: float = 1e-12) -> bool:
    """解析判定：线段 AB 是否与球 (center, r) 相交/相切。"""
    A = np.asarray(A, float); B = np.asarray(B, float); C = np.asarray(center, float)
    d = B - A
    a = float(np.dot(d, d))
    if a <= eps:
        return float(np.linalg.norm(A - C)) <= r + eps
    f = A - C
    b = 2.0 * float(np.dot(f, d))
    c = float(np.dot(f, f) - r * r)
    disc = b*b - 4*a*c
    if disc < -1e-12:  # 明确小于 0
        return False
    if disc < 0.0:
        disc = 0.0
    sqrtD = disc ** 0.5
    u1 = (-b - sqrtD) / (2*a)
    u2 = (-b + sqrtD) / (2*a)
    if 0.0 - 1e-12 <= u1 <= 1.0 + 1e-12 or 0.0 - 1e-12 <= u2 <= 1.0 + 1e-12:
        return True
    # 整段在球内（极少见）：c<=0 表示 A 在球内；还需看 B
    if c <= 0.0:
        return True
    return False